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Abstract 

A cluster algorithm is presented for the Z2 Kalb-Ramond plaquette 
model in four dimensions which dramatically reduces critical slowing. The 
critical exponent z is reduced from z > 2 (standard Metropolis algorithm) 
to z = 0.32 ± 0.06. The Cluster algorithm updates the monopole configu- 
ration known to be responsible for the second order phase transition. 



1 Introduction 



The standard local algorithms for updating monte-carlo simulations become 
very inefficient when applied to systems near a phase transition. As the critical 
point is approached correlation lengths diverge signifying large scale structures on 
the lattice. These structures make simulations much more difficult. In order to 
significantly change the relevant objects, lattice variables with large separations 
have to be changed coherently. 

As the linear size (and therefore correlation lengths at criticality) is increased 
the algorithm has to be applied more times to achieve a statistically independent 
configuration. Thus as the lattice size is increased not only will each update take 
longer but many more updates will be needed. 

Cluster algorithms try to address this problem by identifying the structures 
relevant to the dynamics of the phase transition and using them as the dynamical 
variables. In this way it is possible to construct an algorithm which acts directly 
on the relevant objects while still producing configurations with the correct sta- 
tistical weight. 

Cluster methods were ffist introduced for the Ising model 0, ^ and have since 
been extended to other spin models, in some cases reducing the critical exponent 
to near zero. Attempts to apply the same principles to gauge theories have been 
successful for some models such as the Z2 gauge theory in three dimensions 
and some vertex models However simulations of other gauge groups, like 
SU(3), near criticality remain plagued by Critical slowing down. 

The Z2 plaquette theory [^, 0, being the centre of SU(2) is thought to be 
responsible for the phase transition seen in its adjoint representation, and there- 
fore closely related to the second order point in the fundamental- adjoint phase 
diagram p, H. 

In 0] Ben-Av et.al. showed that the dynamics governing the phase transition 
in a three dimensional Z2 gauge theory is governed by gauge independent vor- 
tex loops and developed an algorithm to stochastically update these loops and 
thus dramatically reduce critical slowing down. Here we will show that these 
techniques can be extended to the 4 dimensional plaquette theory where the 
dynamics are very similar. 



The dynamical variables Px.^iv are defined on the plaquettes of a 4 dimensional 
hyper-cubic lattice and can take on values ±1. The action is defined on 3 dimen- 
sional hypercubes C and is given by: 



2 The Model 



S = ^Se((5C) 



(1) 



c 



where 



Sc(5C) 



(2) 



<5C 
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and 6C are the plaquettes on the surface of the cube C . 

Each 3-cube has a value of ±1 and corresponds to a hnk on the dual lattice. 
Configurations can be uniquely defined, up to a gauge transformation by the links 
on the dual lattice. The bianchi identity for this model constrains the dual lattice 
configuration to be a system of frustrated (-1) dual link loops embedded in a sea 
of satisfied (+1) dual links. These dual loops are the world lines of monopoles 

This conserved monopole current M" is given by 

M; = e'^^^^ApP,,^, (3) 

with 

a^m; = (4) 

Here lies on links of the dual lattice and z labels the dual site corresponding 
to the hypercube at x. Integrating Px,^j.u around a 3-cube gives 

M: = j>Px^^,dC (5) 

Summing this over all 3-cubes (or dual links) gives the monopole world line 
density which is equivalent to the Gibbs free energy. 

This model is dual to the 4 dimensional Ising model which has been extensively 



studied |]T0[ and gives the most accurate value for the Z2 phase transition at 



13 = 0.953. 



3 The Algorithm 

The algorithm used is a single cluster 0] update performed on the dual lattice 
(ie on the 3-cubes) and is therefore gauge invariant. It follows the same structure 
as the 3 dimensional case 0]. 

The algorithm works by modifying the path taken by the monopole world 
lines using the Boltzmann factor to weight any changes to the configuration. A 
random starting point is chosen to build a tree of possible paths. This tree is 
called the graph of deletions. Fig.^ gives an example of how such a tree would 
be constructed (the example has been restricted to 2 dimensions for clarity). 
Starting with the previous configuration pick a random point to be the root 
(Fig.|I|a). Visit each direction and add the neighbouring node with probability 

P, = e'^(^-i) (6) 

Thus if an existing loop is encountered it is automatically added to the graph 
of deletions. Non-frustrated dual links are added with a probability dependent 
on the coupling /?. (Fig]^) shows a possible graph of deletions. All dual hnks 
not in the graph of deletions are frozen. 

The bianchi identity can be used to uniquely update the external legs of 
the tree as valid configurations can only contain closed loops (Fig.|I]c). The only 
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Figure 1: a) The initial monopole configuration and a randomly picked root, b) A 
graph of deletions built from the root, c) The three loops which can be updated 
independently. 
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(a) L1,L2and L3 'On' 



(c) L1,L2'0n', L3'0ff' 



(b) L1 and L3 'On', L2 'Off 
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(d) L1 'On', L2 and L3 'Off 



Figure 2: 4 of the 8 configurations possible from tliis graph of deletions. The 
other 4 have Ll turned off. 



degrees of freedom left for changing the deleted links is a random choice of turning 
each loop on or off. That is choosing the links in each loop to be frustrated (on) 
or satisfied (off). Fig.|^ shows 4 of the new configurations which are possible by 
updating the graph of deletions. The remaining 4 possibilities have Ll turned 
off. Each outcome is weighted equally. 

Using the language of Kandel and Domany |ll[| the dual links forming the 
cluster are said to be deleted and form a graph of deletions. All other dual 
links are frozen. In this language the graph of deletions is updated by randomly 
choosing a new monopole configuration which does not change any of the frozen 
dual links. 

It is easy to see how the algorithm can modify or even destroy an existing 
monopole loop while also being able to create new loops of arbitrary size. 
The algorithm can easily be implemented on a computer as follows. 

1. Pick a random dual site (Hypercube) to be the root of the graph of dele- 
tions. 

2. Visit each dual link (3-cube) connected to this site and add it to the graph of 
deletions with probability Pd = e^^^^'^^ These links are said to be deleted. 

3. Keep repeating the above step for each of the newly added dual sites until 
all of the dual links connected to the graph of deletions have been visited 
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making sure that each dual hnk is only ever visited once. 

4. Identify all the loops within the graph of deletions and choose one dual link 
for each loop to be a free link. This leaves a spanning tree. 

5. Randomly assign a value to each of the free links. 

6. All remaining links in the spanning tree are constrained so the bianchi iden- 
tity can be used to update each link in the tree and create a valid configu- 
ration. 

The algorithm is obviously ergodic as there is a finite probability of deleting 
all of the dual links and choosing a completely random configuration from the 
set of all valid configurations. Detailed balance also follows along the usual lines. 



4 Measurements and error analysis 

To measure how effective the algorithm is at producing independent configura- 
tions we calculated the normalised autocorrelation function for the Monopole 
density M. 

Integrating (|^ gives the integrated autocorrelation time r^nt ® . 

+00 

Tint,M = ^ X! (8) 
^ t=-oo 

Tint measures how many more updates are needed due to the configurations not 
being independent. For a sample mean: 

1 " 

M = -J2Mt (9) 

the variance in M is 

Var(M) = -2Tint(< > - < M >'^) (10) 
n 

for n >> r. This is 2rj„t times larger than expected for independent configura- 
tions. 

The Autocorrelation function decays exponentially at least for large s. Fitting 
the tail of (|^ to an exponential gives Texp the exponential autocorrelation time. 
Tint is defined so as to make Tint ~ ^exp for large r. The autocorrelation time 
diverges as a power of the correlation length so at the critical point r varies with 
the lattice size L as: 



T^L' (11) 

The critical exponent z quantifies the extent to which critical slowing down 
is affecting the simulation. Local algorithms tend to have z 2 for models with 
a phase transition. 
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5 Results 



Simulations were performed on Dec 3400AXP Workstations. For comparison a 
standard Metropolis algorithm was also studied. The monopole density was 
calculated after each update. A number of simulations were performed for /3 
values around the critical point in order to find the largest autocorrelation time 
for the lattice size. Both exponential (Texp) and integrated (Tint) autocorrelation 
time were measured for the Monopole density. Both were found to be in agree- 
ment although for the cluster algorithm Texp proved easier to calculate as the 
autocorrelation functions fitted a exponential very closely before the tail became 
too noisy. Errors in the autocorrelation function were calculated using blocking. 
We noticed auto-correlation functions for metropolis were not pure exponentials. 
To fit this to an exponential decay it was necessary to find a window in the 
autocorrelation function after an initial power law dependence but before noise 
takes over. Metropolis gave a critical exponent of > 2 as expected for a local 
algorithm. 

As our algorithm only updates one cluster at a time the effect of each appli- 
cation can vary from updating no links to updating the entire lattice. In order 
to make a useful comparison of the CPU requirements for each algorithm the 
sweeps of the cluster algorithm were scaled by a factor C representing the aver- 
age cluster size. This factor was calculated as follows. Each application of the 
algorithm is said to be one hit. This hit touches C dual links of the lattice. S 
of these links form the spanning tree. F are free links. R links are looked at but 
rejected. The graph of deletions therefore contains S + F links. Fig.l shows the 
CPU load is roughly proportional to the cluster size C which we have defined as 
C = S + F + R. 

One sweep of the lattice is defined as N/C hits where N is the total number 
of links on the lattice. Thus a sweep takes approximately the same time for all 
(3 and each dual link is touched on average once per sweep. Each Cluster sweep 
took approximately the same time as a metropolis sweep. The extra complexity 
is compensated for by working on the smaller dual lattice where it is easier to 
calculate the action. Our code for metropolis updated approximately 0.42 x 10^ 
plaquettes/s a second while the cluster algorithm updates approximately 0.3 x 10^ 
dual links/s. 

Runs consisted of 1 x 10^ — 2x10^ clusters which gives approximately 6x10^ 
— 6 X 10^ sweeps. Fig ^ shows the results for the autocorrelation time Tint against 
lattice size. The solid lines are our best fits of Zc = 0.32 ± 0.06 for the cluster 
algorithm and Zc = 2.58 ± 0.1 for metropolis. 

6 Conclusions 

The Cluster algorithm presented above becomes more efficient than Metropolis 
for lattice sizes as small as 4^ and is approximately 1000 times more efficient for 
a 15'^ lattice. Comparing our algorithm to Wolff's for the 4D Ising model gives 
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Figure 3: A comparison of Cluster sizes and CPU requirements for a 5^ lattice 
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Figure 4: Monte-carlo simulations of Z2 Kalb-Ramond model in 4D using 
Metropolis and Cluster algorithms. Autocorrelation time vs Lattice size. 
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Model 


Swendsen-Wang 


Single Cluster 


Z2 4D 




0.32 ±0.06 


Ising 4D 


0.83 ±0.01 


0.29 ±0.07 


Z2 3D 


0.73 ±0.06 




Ising 3D 


0.75 ±0.01 


0.39 ±0.01 



Table 1: Dynamical Critical Exponents for Z2 models and the Ising models to 
which they are dual. 



further support to the idea of universality classes. You can see from Table |T] that 
dynamical critical exponents for the models and their corresponding dual Ising 
models are the same within errorbars. 

The embedding of Ising spins has been used successfully to extended the 
realm of cluster algorithms to other spin systems although not all embeddings 
work It is hoped that embedding this Z2 algorithm in a full SU{2) theory 
will give some insight into the structure of the SU (2) fundamental-adjoint phase 
diagram H & O, |T3|, |T3| and the bulk transition of pure SU(2). 
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